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ABSTRACT 

We investigate the effect of black hole spin on warped or misaligned accretion discs — in par¬ 
ticular i) whether or not the inner disc edge aligns with the black hole spin and ii) whether the 
disc can maintain a smooth transition between an aligned inner disc and a misaligned outer 
disc, known as the Bardeen-Petterson effect. We employ high resolution 3D smoothed parti¬ 
cle hydrodynamics simulations of a-discs subject to Lense-Thirring precession, focussing on 
the bending wave regime where the disc viscosity is smaller than the aspect ratio a < H/It. 
We first address the controversy in the literature regarding possible steady-state oscillations of 
the tilt close to the black hole. We successfully recover such oscillations in 3D at both small 
and moderate inclinations (< 15°), provided both Lense-Thirring and Einstein precession are 
present, sufficient resolution is employed, and provided the disc is not so thick so as to simply 
accrete misaligned. Second, we find that discs inclined by more than a few degrees in general 
steepen and break rather than maintain a smooth transition, again in contrast to previous find¬ 
ings, but only once the disc scale height is adequately resolved. Finally, we find that when the 
disc plane is misaligned to the black hole spin by a large angle, the disc ‘tears’ into discrete 
rings which precess effectively independently and cause rapid accretion, consistent with pre¬ 
vious findings in the diffusive regime (a > H/R). Thus misalignment between the disc and 
the spin axis of the black hole provides a robust mechanism for growing black holes quickly, 
regardless of whether the disc is thick or thin. 

Key words: accretion, accretion discs — black hole physics — hydrodynamics — galaxies: 
jets — (galaxies:) quasars: supermassive black holes 


1 INTRODUCTION 


iBardeen & Pettersonl Il975h first computed the evolution of a 
warped accretion disc subjected to Lense-Thirring precession 
jLense & Thirringjl 191 8h caused by frame-dragging from the spin 
of a central black hole. Although their original equa tions were later 
shown to be incorrect JPanaloizou & Pringlei 1 198 31 ) . their qualita¬ 
tive findings of an aligned inner disc smoothly connected to a mis¬ 
aligned outer disc — the ‘Bardeen-Petterson effect’ — has been 
confirmed both from subsequent ID calculation s with corrected 
equations l iKumar & Pringleiri985i : lPringlej|T992h and in three di- 
mensional smoothed particle hydrodynamics (SPH) simulations by 
INelson & Papaloi zoul 1 200d). 

IPapaloizou & Pringle! 11983lf showed that there were two 
regimes for warp propagation in discs depending on the ra¬ 
tio of the effective viscosity a to the aspect ratio H/R. For 
a > H/R, warps can be described with a diffusion equation, 
whereas for a < H/R they propagate as bending waves at half 
the sound speed (iPanaloizou & LiiJ 119950 . Previous studies by 
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iKumar & Pringlei d 19851) and IPringlel (1 19921) were performed in 
the diffusive regime. The first studies performed in the bending 
wave regime J 1 vanov & Illarionov! 1997l : lLubow. Qgilvie & Pringlei 
I2002L hereafter LOP02) suggested a conflict with the Bardeen- 
Petterson picture — finding that the black hole spin could drive the 
disc tilt into a steady state that is oscillatory and non-zero. This im¬ 
plies that the inner edge of the disc may be misaligned with respect 
to the black hole spin. As simple pictures of gas accretion favour 
alignment of the outer accretion disc with the galaxy disc, this has 
been suggested as an explanation for the observed random orienta¬ 
tion of j ets with respect to their host galaxies dKinnev et alj|200ol : 
1LOP02I) . However, unlike the diffusive regime for which there now 
exists a full non-linear the ory de scribing warps of arbitrary ampli¬ 
tude and a dOgilviel lT999. 2000) only linear theory exists for the 
bending wave regime IPanaloizou&Linl l 1 9951 : iLubow & Qgilviel 
l2QQ0i:lLOP02t though see lOgilviel2006h : hence these studies apply 
only to small amplitude warps and could not account for non-linear 
effects. 


The first 3D numerica l simulations of the B a rdeen -Petterson 
effect were performed by INelson & Papaloizoul J2000l ). using a 
post-Newtonian description of the central potential. Their simu¬ 
lations in both the wavelike and diffusive regimes largely con¬ 
firmed the Bardeen-Petterson effect, namely an aligned inner edge. 
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a smooth transition to the outer disc and importantly, no evi dence of 
oscillations in t he tilt. The reason for the disc repancy with lLQP02l 
is unclear, with iNelson & Papaloizoul J 200(1) claiming that ‘non¬ 
linear effects lead_to the damping of these short wavelength fea¬ 
tures’ whilst IlOPQ2 | show that these features are not small com¬ 
pared to the local disc scale height. Further complicating matters, 
recent simulations that included a full general relativistic (GR) 
treatment may have found these oscillations and a non-z ero tilt at 
the inner edge dFrasile et al.ll2007l :lzhura vlev e t al l2014h . 

A further challenge to the Bardeen-Petterson description has 
arisen from recent simulations of warped discs in the diffusive 
regime (a > H/R). 3D simulations of warps in isolated discs 
have demonst rated close agreement with the non-linear theory of 
lOgilviej j 1999h (see lLodato & Pricel2010l) . However, in simulations 
of warps driven by Lense-Thirring precession it was found that in¬ 
stead of maintaining a smooth transition (as described by the linear 
Bardeen-Petterson effect), discs break when the angle b etwee n the 
disc a nd the black hole spin is more than a few degrees dNixon et al .1 
1 201 2bl) . Rings of gas were found to be tom off the disc, process¬ 
ing effectively independently before being accreted. The net effect 
was a much higher accretion rate on to the black hole. Thus the 
Bardeen-Petterson idea of a smooth transition does not appear to 
hold in the diffusive regime for an arbitrary choice of parameters. It 
is not clear whether or not misaligned discs in the wavelike regime 
will also break. 


Here we reexamine the Bardeen-Petterson effect in black hole 
accretion discs. Building on the success of earlier stud ies we u se th e 
SPH c ode PHANTOM t o model warped discs i n 3D dLodato & Price 


201(1 

Nixon 20121: Nixon. Kins & Price 201 2al; Nixon et al. 

2012 H; 

Facchini, Lodato & Pric3 2013; Nixon, Kins & Pricel2013l: 

Martin et all 2014a. Martin et all 2014H). Here we focus on the 


wavelike propagation regime (a < H/R). We start by consi der¬ 
ing the possible reasons for the discrepancy between IlOPQ 2 | and 
[Nelson & Papaloizoul (200y), as well as two possible complications 
to the Bardeen-Petterson picture of an aligned inner disc smoothly 
joined to a misaligned outer disc in Section 0 We present the nu¬ 
merical method and tests of wavelike warp propagation with SPH 
in Section[3] In Section[4]we present a suite of 3D simulations de¬ 
signed to confirm under what circumstances the Bardeen-Petterson 
description holds. In Section[5]we discuss the implications of these 
and in Section[ 6 ]we outline our conclusions. 

An obvious caveat is that we do not consider magnetic 
fields, even though it is widely accepted that magnetorota- 
tional instability ( MRI) is the contro lling mechanism for vis¬ 
cosity in the disc dBalbus & Hawlevl ll99ll) . However, we can 
still capture the dominant behaviour in the disc as magnetic 
fields have be en shown to h ave little e ffec t on the geometri¬ 
cal evolution dSora thia. Kr olik & Hawlevl 120131) . We also follow 
INelson & Papaloizoul d2QQol) ' in using using a post-Newtonian ap¬ 
proach instead of general relativity (GR). As we will see, one of the 
findings of this paper is that this approximation must be considered 
carefully in order to capture the combination of relativistic effects 
that lead to tilt oscillations. 


2 DOES THE BARDEEN-PETTERSON EFFECT HOLD 
IN THE WAVELIKE REGIME? 

We consider three possible ways that the Bardeen-Petterson effect 
may be violated in wavelike discs. Firstly, radial oscillations in the 
tilt of the disc may prevent the disc from aligning at the inner edge. 
Secondly, the smooth transition between aligned and misaligned 


material may be broken if the disc tears, as has been observ ed in the 
diffusive regime dNixon & Kind2012l : Nixon et alj2012bi) . Finally, 
it may not be possible for the disc to find a steady state if the disc 
is relatively thick and the viscous time is short. 


2.1 Is the inner disc aligned? 

1LOP02I considered warps in geometrically thin, almost Keplerian 
discs described by a surface density E(-R) and angular velocity 

Q. (R). The scale-height of the disc is given by H(R) = c a /H, 
where Cs(R) is the sound speed in the disc. Their description is 
one dimensional in the sense that the total angular momentum in 
the disc L is a function only of the cylindrical radial coordinate, 

R. The disc is then discretised into a series of rings, each de¬ 
scribed by the orientation of its tilt and twist angle. The tilt an¬ 
gle /3 is measured from the 2 -axis, and if this angle varies with 
radius the disc is considered to be warped. The twist angle 7 is 
measured from an axis that is perpendicular to the 2 -axis, and sim¬ 
ilarly, if the twist angle varies with radius the disc is twisted. These 
two angles can be related to the unit angular momentum v ector by 
l = L/L = (cos 7 sin /3, sin 7 sin /3, cos /3) (Prin rfil996l) . 

We assume an a disc viscosity where v = ac s H 
dShakura & Sunvaevll 1 9731) . For accretion d iscs with a < H/R, 
the warp propagates as a disp ersive wave JPanaloizou & Pringle! 
1 19831 : Papaloizo u & Lin| 1 1 995D . Assuming that the disc is nearly 
Keplerian and not self-gravit ating, the equations of motion describ- 
ing the wave propagation are dLubow & Qgilviel200d : lLubow et ahl 
1 20021) 




( 1 ) 


dG 
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( 2 ) 


Here k is the epicyclic frequency, G represents the internal 
horizon tal tor que in the disc and T is the external torque per unit 
area. 1 lopo 2 < chose a complex representation where the warp is 
given by W = l x + il y and the internal torque as G = G x + iG y . 
This allows Equations Q] and [2] to be rewritten as 
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These equations describe the propagation of a warp in the linear 
regime, and were solved numerically by iLOPOl to find the steady 
state shape of the disc around a Kerr black hole. In this case the 
apsidal and nodal precession frequencies in the disc (scaled by fl) 
can be approximated to first order from the Kerr metric as dkatol 
11990h 


V LOP = 


k 2 - n 2 

2 H 2 


2 R ’ 


(5) 


. n 2 -n 2 a fR s \ 3/2 

ClOP “ 2D 2 _ y/2 V R ) 


(6) 


where R s = 2GM/c 2 and a is the black hole spin. These fre¬ 
quencies are used in the solution by inserting the m dire ctly into 
Equations[3]and[4] An example of the solution bv lLOPQ^ is shown 
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Figure 1 . The steady state tilt profile found by |LQP02l as in their Figure 3 and shown to the same time (at eight equally spaced times). The original apsidal 
and nodal frequencies are used in the left panel (Equations |5]and|6) with the same sign. In the right panel we have reversed the sign of the nodal frequency so 
the precession frequencies have different signs, and there are no oscillations present in this steady state (similar to lLOPojl Figure 6). Here the tilt is shown as 
a fraction of the maximum tilt, the initial condition is shown with the dashed line and fii n = 4 R g . 


in the left of FigureQ] with the same parameters used in their work. 
The steady state solution is formed from the interaction of the ingo¬ 
ing and the outgoing bending waves, where the outgoing waves are 
created b y the reflection of the ingoing waves at the inner boundary 
1LQP02 ). Here the oscillatory behaviour of the steady state near the 
inner edge is clear, as is the non-zero tilt at the inner edge. 

It is known that the relative signs of the apsidal and nodal 
frequencies determines wh ether the solution is oscil latory o r not 
dlvanov & Illarionov! 19971) . The frequencies used bv lLOPO^ have 
the same sign, leading to radial oscillations in the steady state tilt 
profile. We confirm that the oscillatory profile is dependent only on 
the signs of the frequencies by changing the s ign of C lop (equiva¬ 
lent to modelling a retrograde black hole, see lLOP02L Figure 6) in 
the right hand panel of Figure Q] The two solutions evolve in the 
same manner with the exception of the oscillations near the inner 
edge. 

While one is free to set the precession frequencies directly 
when solving Equations [3] and [4] in 3D the nodal precession can 
be induced directly (e.g. using the post-Newtonian description of 
Lense-Thirring precession from a spinning black hole) and the ap¬ 
sidal precession (i.e. Einstein precession) arises indirectly from the 
central potential. Hence, it is possible for the choice of poten¬ 
tial to preclude oscillations from the steady state solution in 3D 
simulations. It is then not surprising that simulations that do not 
take the apsidal pre cessi on int o accoun t as above also do not re- 
port tilt oscillation s JSorat hia et alj2013l) . However, simulations by 
iNelson & Papaloizoul d20oT did make use of a potential that re¬ 
sulted in apsidal and nodal precession frequencies with the same 
sign but did not find oscillations. Here we use high resolution sim¬ 
ulations along with a potential that leads to precession frequencies 
of the same sign to investigate this discrepancy. 


2.2 When does the disc break? 

The derivation of EquationsQ]and[2]assumes that the inclination of 
the disc is linear. From previous results in the non-linear regime we 
would anticipate that the disc may break when the external torque 
applied to the disc is stronger than the internal torque. Here the 
internal disc communication is governed by a combination of pres¬ 


sure and viscosity. The viscous torqu e that acts between successive , 
discrete rings in the disc is given by dLvnden-Bell & Pringlejl 19741) 

G = 3ttzz£ (GMR) 1/2 . (7) 


Lense-Thirring precession causes the rings that make up the disc 
to precess. Per uni t area on the disc, this torque is given by (e.g. 
iNixon et ~akll2012bl) 

T = ^” SjR2f2 l Sin6l l (if) 5 (8) 

where R g = GM/c 2 , a is the black hole spin and 6 is the angle 
between the plane of the disc and the direction of the black hole 
spin. If the external torque applied to the disc is greater than the 
internal torque maintaining the disc, the rings will precess inde¬ 
pendently _faster than the disc is able to communicate the preces¬ 
sion l INixon et al .1120120 ). This will result in the disc being sep¬ 
arated and breaking, perhaps into differentially precessing rings. 
Assuming that the disc has no initial warp and that internal com¬ 
munication is dominated by viscosity, a comparison of the above 
torque s predicts a maximum radius that it is possible for this to 
occur (Ni xon e t al J l2013ll 


Rbreak 


4 a 
3 a' 


2/3 


R r . 


(9) 


This approximate relationship places an upper bound on the break¬ 
ing radius of the disc at a given angle. For the typical parame¬ 
ters used in this paper, we have H/R = 0.05, R ou t = 407?i n , 
a = 0.01 and a = 0.9. At the outer edge of the disc, the above 
relation then reduces to 

7?break < 45f? in (sin d) 2/3 R g . (10) 


This predicts that tearing may occur in the disc for inclina¬ 
tions of more than 6°. At this inclination or greater one would ex¬ 
pect the discs to break rather than align. However, in the bending 
wave regime that we consider here, the internal communication is 
dominated by pressure. In this case we can estimate the radius at 
which the disc will break by comparing the sound crossi n g and 
the precession timescales in the disc. Following lNix o n et al. 120131) 
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and assuming that the disc is inviscid (and hence not taking into 
account any wave damping) we find that 

-Rbreak.t < ( 4a| sin 6 »| —J i? g . (11) 

We therefore only expect the disc to break closer to the black hole 
than Equation 1101 


2.3 Can the disc accrete misaligned? 


A further assumption made in developing Equations Q] and [2] was 
that the viscous timescale in the disc is much larger than any 
other timescale, equivalent to assuming that the disc is replen- 
ished from ra dii outside the computational domain, or a <C H/R 
dLubow et al.ll2Qo3) . This implies that the surface density profile 
does not change during the evolution of these equations, which 
is valid until the warp reaches the outer boundary. We can quan¬ 
tify this approximation during the evolution of the equations by 
considering the ratio of the wave and visc ous timescales (as in 
iLodato & PringleflilOOdlFacc hini et ~aPl2013h 


twELVe 2 Ru n H 
t v ~ c s R 2 " R ' 


( 12 ) 


For a < H/R <C 1, the above relation implies that the viscous 
time is much greater than the sound crossing time that t he warp 
communicates, so we can neglect mass accretion. Indeed, IlOPCQ 
neglected the evolution of the surface density profile completely in 
their solution, equivalent to assuming no mass is accreted at all. 
However for their disc H/R = 0.1030, and so it is not clear that 
this assumption holds. Additionally, the viscous time can be written 
as 


=-(Ay 2 

afi \ RJ 


(13) 


In this form, it is clear that increasing the aspect ratio of a disc re¬ 
sults in a significant decrease in the viscous time. At a given radius 
R, when the viscous timescale is comparable to (or smaller than) 
the precession timescale at that same radius, accretion dominates. 
In thick discs (or tori) it may then be possible for the material in 
the outer disc to be accreted before it has a chance to align. The tilt 
profile in this case will not reach a steady state but instead be deter¬ 
mined by the inward flux of angular momentum. Thus in relatively 
thick discs no tilt oscillations would be expected. 


3 NUMERICAL METHOD 

We use 3D simulations to investigate whether the Bardeen- 
Petterson effect holds in accretion discs focussing on tilt oscilla¬ 
tions, large inclinations and misaligned accretion in the limit where 
a < H/R. We use PHANTOM, an efficient and low memory SPH 
code iLodato & Price! [2oTcl : IPrice & Federratl]|20 1 (1 |Pricell201~3) . 
This code has been widely used to model warped discs in the diffu¬ 
sive regime and to investigate related problems with warped wave¬ 
like discs (see Section[Tjl. The physical disc viscosity in the disc is 
represente d in the code using the artificial viscosity (cyav) method 
outlined in lLodato & Price! J20 1 (ill . 


3.1 Modelling Lense-Thirring Precession 

The Lense-Thirring precession exerted by the black hole is 
modelled using a post-Newtonian approximation, represented 


as a first order (in v/c) correction in the momentum equa- 
tion. Taking this into accou nt, the momentum equation becomes 
( iNelson & PapaloizoulbOQOh 

37 = —i\7P + v x A -V$ + S’vtsc. (14) 

dt p 

Here S v isc is the viscous force per unit mass and v x h is the gravo- 
magnetic force per unit mass, with h given by 


h = 2S _ 6(5 -r)r 
R 3 R 5 


(15) 


and 5 = aG 2 M 2 k/c 3 , pointing in same direction as the black hole 
spin. 

Assuming a thin disc with linear disturbances, the quantity 
v x h can be expressed as 


v x h = 



6 S(zr x v) 

+ IP 


(16) 


This expression is mainly useful for setting up the initial orbital 
velocities for the disc particles (Section[T2]l- 


3.1.1 Implementation in Leapfrog 

A complication to the implementation in the code is that we use 
a leapfrog integrator in the ‘Velocity-Verlet’ form, where the po¬ 
sitions and velocities of the particles are updated from time t n to 
t n+1 according to 

v "+5 = v n + ^Ata n , (17) 

x n+1 = x n + Atv n+ ?, (18) 

v n+1 =v n+ i + i Aia n+1 . (19) 

However, the acceleration caused by Lense-Thirring precession de¬ 
pends on velocity. Thus CD becomes implicit. This can be easily 
solved by writing the corrector step in the form 

v" +1 = v n+ i + t (v n+1 x h n+1 ) , (20) 

where fip ^ 1 contains the position-dependent terms. This forms a 
set of three linear equations for each component of v n+1 , that we 
solve analytically by inverting the resulting 3x3 matrix. 


3.1.2 Potentials 


We make use of two gravitatio nal potentials in this work. The 
first was previously introduced bv lNelson & Papaloizoul ( l200d) . re¬ 
ferred to as the Einstein potential (see their Equation 8 ). In our no¬ 
tation it is given as 


$e (R) 


GM / 3-R g \ 

r ( + R ) ' 


( 21 ) 


This potential was introduced because it prevents the gravitational 
force tending to infinity as the radius decreases. However, it also re¬ 
sults in the correct apsidal precession frequency at large distances 
from the black hole and has the same sign as the nodal frequency 
( INelson & Papaloizoul200d) . This is in contrast to the standard Ke- 
plerian potential 


HR) 


GM 

~R~' 


( 22 ) 


The standard potential 1221 was used in all of the non-linear simu¬ 
lations, except for Figures[4}{3 w liere (|2 1 1 was used. 
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3.1.3 Precession Frequencies 

We calculate the apsidal and nodal precession frequencies in our 
disc using the standard (Newtonian) definitions for the epicyclic 
and vertical frequencies, 

k 2 = 4 n 2 +R^-(Q 2 ), (23) 


n 


2 

z 


d 2 4>(R) 

dz 2 


2 = 0 


(24) 


Following [Nelson & Papaloizoul (|2000|) . we compute these using 
an effective potential that takes in to account the second and third 
terms of Equationll4l 


4> e s{R) 


HR) + 


4 sVGM 

5f ? 5 / 2 


6 Sz 2 VGM 

R 9/2 


(25) 


This potential accounts for both the normal Keplerian potential and 
the correction due to the v x h term, where &(R) could be repre¬ 
sented by either Equation[2T]or[22] Firstly considering the Einstein 
potential, using Equations [23] and [24] the post-Newtonian apsidal 
and nodal precession frequencies (scaled by fl) are given by 


2 fl 2 


6 GMR g 
R* 


ssVgm' 

R 9/2 


(26) 


Simulation 

on 

a 

a 

OiAV 

PSl 

30 

0.1 

0.01 

0.395 

PS2 

30 

0.1 

0.03 

1.186 

PS 3 

30 

0.3 

0.01 

0.395 

PS4 

30 

0.3 

0.03 

1.186 

PS5 

30 

0.5 

0.01 

0.395 

PS6 

30 

0.5 

0.03 

1.186 

PS7 

30 

0.7 

0.01 

0.395 

PS 8 

30 

0.7 

0.03 

1.186 

PS9 

30 

0.9 

0.01 

0.395 

PS 10 

30 

0.9 

0.03 

0.186 

A1 

0 

0.9 

0.01 

0.395 

A2 

15 

0.9 

0.01 

0.395 

A3 

30 

0.9 

0.01 

0.395 

A4 

45 

0.9 

0.01 

0.395 

A5 

60 

0.9 

0.01 

0.395 

A6 

90 

0.9 

0.01 

0.395 

A7 

120 

0.9 

0.01 

0.395 

A8 

150 

0.9 

0.01 

0.395 


Table 1. Simula ti on p arameters, including the spin (a) and 
IShakura & Sunvaevl 41973T) viscosity (a) and artificial viscosity ( olav )• 
Unless otherwise noted, the accretion discs also had H/R = 0.05, an 
outer radius of 40-Rin and made use of 10 7 particles. 


_ -2 SVGM 

where Cl 2 is given by 

o2 GM ( 67? g \ 2 SVGM 

Qe -~r ^{ 1 + ^~)~ RV 2 


(27) 


(28) 


As the signs of the apsidal and nodal precession frequencies here 
are the same throughout the disc, this potential will allow an 
oscillatory profile to develop. Considering now a standard post- 
Newtonian potential, given by Equation[22] we find the apsidal and 
nodal precession frequencies to be (again scaled by fl) 


?7pn 



(29) 


-45 

Cpn = — ,_ _. (30) 

2VGMR 3 / 2 - 45 

Here w e note an important difference with respect to the solution 
used bv lLOPO^ . As the signs of the precession frequencies here are 
opposite, the steady state tilt profile will not have oscillations if the 
potential in Equation[22]is used. 


3.2 Initial Conditions and Scope 

Unless otherwise stated, the discs presented all made use of 10' 
particles, and simulations with 10 6 and 10 5 were also conducted 
to check convergence. We note that each time the resolution is 
changed between simulations with otherwise the same parameters, 
the artifici al viscosity is a ltered according to the scaling described 
in iLodato & Pricel QQloh so that the discs have the same a inde¬ 
pendent of which resolution is used. The locally isothermal sound 
speed in the disc was set to Cs(R) = c s ,m(-R/-Ri n ) -9 and the sur¬ 
face density profile E(J?) = Ei n (i?/i?i n ) _p , where p = 3/2 and 
q = 3/4 to give a constant a vis cosity in the disc and uniform 
resolution dLodato & Pringleil2007r) . Each disc was initially set up 


aligned to the black hole spin, with the particles arranged using a 
Monte Carlo placement method. Each particle was then rotated by 
the inclination angle and assigned a velocity according to the fol¬ 
lowing expression derived from Equationll6l 




R 3 

+ ^r a 


cos($), 


(31) 


where «k is the Keplerian orbital velocity. The discs were there¬ 
fore initially tilted to the black hole spin, but not warped. The re¬ 
sults presented below have time shown in orbits at the inner edge 
and show the tilt as a function of radius only. This was found 
from the simula t ions u sing the method outlined in Section 3.2.6 of 
ILodato & Pricel ( bold) where we used N = 300 spherical shells. 
For all of the simulations the in ner radius was set as Ri n = 4R g , 
in order to compare to the lLQPoj ID code. At small radii we note 
that the absence of GR limits the validity of our results, and indeed 
the need to carefully account for relativistic effects is one of our 
findings. 


3.3 Test of wavelike warp propagation 

To confirm that we can correctly describe the propagation of 
warps in the wa v elike regime, we use the test described by 
[Fragner & Nelson! d20ld) . They simulated a wavelike accretion 
disc with a point mass potential and compared to a ID calcula¬ 
tion, finding agreement at the ~ 10% level. We choose to compare 
to this ID solution instead of the solu ti on fr om Equations Q] and [2] 
because the linear solution from iFragner & Nelson! d2010[) allows 
the surface density to evolve, as occurs in our 3D simulations. 

We conduct a simulation using the same parameters as cited 
in their Figure 1, with H/R = 0.03 and a = 0.001 and an initial 
disturbance of 5°. In this case we do not drive the evolution with 
Lense-Thirring precession, so that we can isolate the behaviour due 
to warp propagation only. Our results are shown in Figure [2] using 
10 6 and 10' particles. As the disc evolves, the disturbance splits 
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Figure 2. Evolution of bending waves in a disc, not subject to Lense- 
Thirring precession. The green (10 6 ) and red (10 7 ) lines show the results 
fr om our 3D si mulatio n. The bl ack lin e shows the results from the ID code 
of iFragner & Nelsoi] (Figure 1 l2010h . using the same initial parameters. 
The agreement between these two solutions confirms that SPH can be used 
to describe the evolution of warp propagation in the wavelike regime. 


into two waves travellin g at half the sound speed (as predicted by 
IPapaloizou & Linl 1995b : one inward and the other outward. By the 
end of the simulation these have fully separated and are beginning 
to interact with the boundaries. 

The SPH solution shows the same behaviour as the ID solu¬ 
tion, and increasing the resolution reduces the discrepancy. How¬ 
ever, at late times the inner edge of the disc there is increasing 
disagreement, most likely due to differences in the inner bound¬ 
ary condition. This test confirms that PHANTOM can be used to 
describe the propagation of warps in the wavelike regime. 


3.4 Test of Lense-Thirring precession 

We also perform a simple test of the Lense-Thirring precession. 
We simulate a disc consisting of test particles with no viscosity and 
zero sound speed (i.e. a = c a = 0) subject to Lense-Thirring pre¬ 
cession. The initial velocities are set using Equation[3T]and the disc 
is inclined at 30°. We then calculate the precession in the disc as a 
function of the radius using the procedure outlined in Appendix lAl 
Figure [3] shows the comparison between the precession mea¬ 
sured from our disc and the predicted precession in the disc, given 
by t p = R 3 /(2a). We find agreement to within measurement un¬ 
certainties throughout the disc. 



Figure 3. Precession timescale measured from an inviscid and pressureless 
3D disc as a function of radius (black circles), compared to the expected 
Lense-Thirring precession (red line). 71; n = 4 R g . 


4 RESULTS 
4.1 Tilt Oscillations 

We first investig ated whether or not the tilt oscillations predicted 
bv ILubow et al. J2002I ) are physical using 3D simulations. The disc 
is initiated with a constant misalignment of 3°, within the linear 
regime required by Equations Q] and [2] We c hose p arameters for 
our simulation similar to that of iLubow et al.1 I I2002I) . with the ex¬ 
ception of the surface density profile, the black hole spin and the 
disc thickness. Additionally, we made use of the Einstein potential 
outlined in Equation[2T] in order to give precession frequencies of 
the same sign (Equations [26] and [27}. We setp = 1.5 and q = 0.75 
so that the disc is uniformly resolved, as discussed in Section [3721 
The disadvantage is that this results in lower amplitude oscillations 
in the ID code. To combat this we encourage larger amplitude os¬ 
cillations by increasing the spin to a = 0.9 and decreasing the disc 
thickness to H/R = 0.05. The evolution of the tilt as a function 
of radius is shown in Figure [4] from a simulation employing 10' 
particles. 

One of t he main differences between t his simulation and those 
conducted by iNelson & Papaloizoul d2000h is the angle of inclina¬ 
tion. While our 3° initial tilt was well in the linear regime required 
by the analytic description in Equation s[j~|and[2l the minimum incli¬ 
nation used by iNelson & Papaloizoul l200(t) was 10°. We explore 
the effect of non-linear inclinations in this potential by misaligning 
the same disc at 15°. Figure [5] shows a cross section of density in 
the inner disc from this calculation. The tilt profile after ~ 600 or¬ 
bits is qualitatively similar to Figure[4] showing the same evolution. 
The quantitative tilt evolution is resolution-dependent, but never¬ 
theless a non-zero tilt and oscillations were found at both medium 
(10® particles) and high (10 7 particles) resolution. 

Figure [6] shows a resolution study of the tilt and surface den¬ 
sity profiles using 10 5 , 10® and 10 7 particles. The main artefact of 
low resolution is that accretion occurs faster and as a result there 
is less mass at the inner edge in the lower resolution calculations. 
Comparison of the surface density profiles indicates that E(iZ) is 
not fully converged near the inner edge, which has a dramatic effect 
on the tilt profiles. However, in all discs there is a non-zero tilt at 
the inner edge and in the 10® and 10' discs radial oscillations are 
observed. The wavelength of these oscillations is consistent with 
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Figure 4. Time evolution o f the a ng le between the disc plane and the black hole spin as a function of radius in a 3D disc subject to Lense-Thirring precession, 
using similar parameters to lLubow et al.1120021) and the same times as in Figure [I] The shape of this profile depends sensitively on the surface density at the 
inner edge and hence on resolution (see Figure[6}. The right panel shows a zoom-in of the inner disc. 



Figure 5. Cross-section view of the steady-state tilt oscillation formed in a disc initially inclined at 15° to the black hole spin (spin axis is vertical with respect 
to the page, i.e. along the z axis). The colour scale shows density. Disc parameters are the same as in Figure]?] but with a larger initial inclination. 


the criteria given by iLubow et al.l d2002h - Even using the preces¬ 
sion frequencies and surface density profiles in the ID code that 
are appropriate to the 3D simulations (Equations 1261 and 127b still 
does not provide a close match with the 3D results. It is not clear if 
thi s discrepancy is due to non-l inear fluid effects, e.g. as discussed 
by iNelson & Papaloizoul d2000l) . or simply requires higher resolu¬ 
tion calculations to obtain numerically converged results. However 
it is clear that with the appropriate potential and system param- 
eters, the disc can display rad ial tilt o scillations as predicted by 
[ivanov & Illarionov! d 19971 ) an d lLOPOl 


Despite the resolution-dependence of our results, we were 
still able to observe tilt oscillations at resolutions used by 
iNelson & Papaloizoul ( hood ) so long as Einstein precession was ac¬ 
counted for. We further investigated whether this might be due to 
the differences in the artificial viscosity parameters used, as we set 
the Von Ne umann-Richtmyer viscosity coef ficient /3av = 2.0 (see 
[Price] 2012h for all of our simulations whilst [Nelson & Papaloizoul 
d200C ) used /3 av = 0. A n onzero /3av visco sity is required to 
prevent particle penetration dMonaghanl 119891) and the absence 
of bul k viscosity is known to be problematic in disc simulations 
dLodato & Pricel 2010h. Thu s with /3av = 0, the simulations of 
INelson & Papaloizoul feood) might not have captured the wave in¬ 
teractions that create the tilt oscillations and the absence of bulk 


viscosity. To check this we co nduct a low-resolution simul ation 
equivalent to simulation El of INelson & Papaloizoul ( 2000 ) and 
/3av = 0. FigureQshows the results (black solid line), compared to 
an equivalent simulation with /3av = 2 (red dashed line) and also 
compared to a higher resolution simulations. At high resolution we 
find tilt oscillations regardless of the value of /3av (green solid and 
blue dotted lines) but we find that using /3av = 0 can indeed erase 
the tilt oscillations at low resolution. The lower panel of Figure [7] 
shows that this is not simply due to the effect on S (R), since two 
of the calculations show very similar surface density profiles but 
rather different evolutions of the inner disc tilt. 


4.2 When does the disc break? 

4.2.1 Bardeen-Petterson Alignment 

A second possible violation of the Bardeen-Petterson picture may 
be that the disc breaks instead of maintaining a smooth transition 
between an aligned inner disc and a misaligned outer disc. In or¬ 
der to investigate this, we simulate a range of discs at 30° whilst 
varying a and a according to the list PS 1-10 in Table 13.01 Here, 
for simplicity, we make use of a standard potential given by Equa¬ 
tion]^] and hence do not expect any oscillatory behaviour. 

Figure [8] shows a 3D rendering of density in one such simu- 
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Figure 6 . Resolution study showing the inclination (tilt angle) and surface 
density as a function of radius at the final time of the disc shown in Figure[4] 
using 10 5 (long dashed green), 10® (short dashed red) and 10' (solid black) 
particles. Increasing the resolution better resolves the surface density profile 
at the inner edge, which strongly affects the final tilt profile found. 

lation with a = 0.9 and a = 0.03 at three different resolutions 
(PS 10). Except for the p otential used, this disc has si milar param¬ 
eters to simulation E3 of lNelson & Papaloizoul l200(j|) which made 
use of 52, 000 particles across a larger radial extent than our sim¬ 
ulations, representing a lower resolution than any of those shown 
in Figure [ 8 ] At our lowest resolution (left panel of Figure [ 8 j, we 
also observe the inner disc aligning and smoothly transitioning to 
an outer, misaligned disc (see their Figure 12). However, at higher 
resolutions this behaviour is no longer observed — the disc instead 
breaks into two distinct sections, with the inner disc aligned with 
the black hole spin and the outer disc remaining misaligned. 

Figure [9] shows that resolving disc breaking is mainly a ques¬ 
tion of resolving the disc scale height. For the lowest resolution 
simulations the resolution length is greater than the scale height of 
the disc and hence disc breaking (on a length scale smaller than H) 
cannot be resolved and a smooth transition is observed. By con¬ 
trast, the two higher resolutions are able to resolve disc breaking. 
Figure [TO] shows the same resolution study performed in Figure [8] 
for all of our discs at 30°, where the green line shows simulations 
that made use of 10 s particles, red shows 10 6 and black shows 10 ' 
particles. The discs are shown after 1500 orbits at the inner edge, 
allowing the warp to propagate all the way to the outer radius. In¬ 
creasing the spin of the black hole increases the rate at which the 
innermost part of the disc aligns. 

Across all of the parameters chosen here, increasing the reso- 
lut ion changes the behav i our fro m the smooth tilt profile observed 
by iNelson & Papaloizoul ( boQOi) to a steepening of the tilt profile 
and ultimately a disc that is broken into distinct sections. The higher 
resolution results show an aligned inner edge, a misaligned outer 
edge and a sharp tilt profile connecting these, representing a break 
in the disc. The discs simulated with lower spins appear to steepen 
and tear faster than those with higher spin as the break occurs fur¬ 
ther out (and hence a longer precession time). This is observed most 
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Figure 7. The effect of bulk viscosity on the final tilt of the accretion 
disc. The solid, black line u ses the same parameters as simulation El of 
iNelson & PapaloizoJ 1200(1) . with 5.2 X 10 4 particles and has no bulk vis¬ 
cosity. Bulk viscosity is used with both the red (short dashes, 5.2 X 10 5 
particles) and green (long dashes, 5.2 X 10® particles) lines. At higher res¬ 
olution and with bulk viscosity tilt oscillations are resolved, but the inner¬ 
most parts of the disc remain unconverged. Here R g = 0.04 and the disc 
has been evolved until the warp has reached the outer edge. 

clearly between the low viscosity, high spin cases. At a = 0.5, the 
tilt steepened and the disc tore before the end of the simulation. For 
the disc with a = 0.7, the tilt began steepening near the end of the 
simulation but was not able to separate, whilst at a = 0.9 the disc 
has not yet begun steepening. We have confirmed that this is the 
case by extending the high resolution simulations of the a = 0.9 
case, and indeed observed steepening to occur at later times. 

As with the previous simulations, FigureQT^demonstrates that 
the simulations are not fully converged, especially when consid¬ 
ering the low spin cases (a < 0.5). For these discs, the discrep¬ 
ancy in the inner tilt is again due to the mass accreted at the inner 
edge of the disc. At low resolutions the inner part of the disc is ac¬ 
creted faster, resulting in less mass near the inner edge. The same 
Lense-Thirring torque then acts on less mass, and is thus not able 
to align the disc to the same extent. At increasing spins this effect 
is observed less, as the higher spin provides a larger torque and so 
even the lowest resolution discs are able to align. In the discs with 
a < 0.5, increasing the resolution leads to a more distinct tear in 
the disc suggesting that our results are consistent. Hence we can 
be confident that these discs do tear, and present an upper limit on 
the radius at which this occurs. As the tearing occurs outside of the 
radius where oscillations were found in Section |4~T1 using similar 
parameters, this behaviour should not be affected by our choice of 
potential. 


4.2.2 Disc Tearing 

To investigate the dependence of disc tearing on the misalignment 
between the disc and the black hole spin in the wavelike regime 
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Figure 8. Structure of the disc with a = 0.9 and a = 0.03 at increasing numerical resolution (left to right). At low resolution the Bardeen-Petterson Effect is 
observed, similar to the results of lNelson & Papaloizoul |200(1 . but at high resolution the disc distinctly tears into two separate sections. The colour indicates 
density, with white being highest. 


we simulated a suite of discs at different inclinations. We again 
make use of the traditional post-Newtonian approximation given 
by Equation 1251 We held a = 0.01 and a = 0.9 constant and 
varied the inclination of the disc between 0° (aligned) and 150°, 
noted in Table l3.13l with Al- 8 . Figure |TT] shows these simulations 
after more than 1500 orbits measured at the inner edge. Each disc 
was initially tilted but not warped. As the simulation progressed, a 
warp evolved in response to the Lense-Thirring torque and in the 
higher inclination cases resulted in the disc breaking. 

At 15° (top right of Figure DU the disc was observed to 
smoothly align to the spin of the black hole. At the end of the simu¬ 
lation, the tilt of the disc was consistent with the Bardeen-Petterson 
effect and is similar to results seen in previous simulations at 10 ° by 
iNelson & Papaloizoul d2000i) . Extending the lower resolution ver¬ 
sion of this simulation (with 10 6 particles) for twice as long shows 
that the disc continues to align with the black hole spin, implying 
that the steady state for this disc is full alignment. Inclining the disc 
at 30° also did not yet result in disc tearing, however this is because 
for this particular choice of viscosity and spin this simulation has 
not been run long enough (see Section l4.2.1l i 

For discs at higher inclinations (> 45°; second, third and 
fourth rows of Figurell It. the inner section of the disc was found to 
align within 50 orbits and a smooth transition was formed between 
this and the outer region of the disc. This transition then steepened 
until the disc broke into two sections that were connected by pre- 
cessing rings of material. Multiple rings of material were torn off 
from the outer, misaligned disc and each was observed to precess 
effectively independently. Towards the end of the simulations, up to 
two rings were precessing at the same time (for example, 120 ° disc 
of Figure mi l and were present for up to ~ 400 orbits. Eventually 
each of these rings settled with and increased the inner, aligned re¬ 
gion of the disc. The disc inclined at 90° (right hand panel in third 
row of Figure [Ill also developed precessing rings of material that 
were accreted. However, for this inclination no inner aligned disc 
was observed and the rings of material were accreted directly onto 
the black hole. 

Figure [T^ shows the instantaneous mass accretion rate by the 
discs at different angles. It can be seen that inclining the disc to the 
spin of the black hole increases the rate of accretion by more than 
an order of magni tude w hen compared to an aligned disc, similar to 
previous findings jNixon et alj2012bh . The discs that form an inner 
aligned disc and precessing rings have even higher accretion rates, 
as the inner disc is continually fed by the rings as they align . When 
taken in context with the results in the diffusive regime jNixon et al.1 
l 2012 bl) . this implies that regardless of whether the disc is thin or 
thick, mass accretion is faster when the disc is inclined. 



Figure 9. Resolution length as a fraction of the disc scale height ((h)/H) 
for the three resolutions (10 s in green long dashes, 10® in short red dashes 
and 10 7 in solid black) shown in Figure[8] As disc breaking occurs on length 
scales smaller than the scale height of the disc, the lowest resolution simu¬ 
lations here cannot resolve breaking behaviour (as the resolution length is 
greater than the scale height throughout the disc). The two higher resolution 
simulations are able to resolve this behaviour. 

Disc tearing has also been observed in the wavelike 
disc r egime for circumbinary discs inclined at high angles 
dFacchini et al.1120131) . In a simulation of a circumbinary disc in¬ 
clined at 60°, their disc separates into two sections and the inner 
one precessed effectively independently of the outer disc. As their 
disc is thicker than ours ( H/R = 0.1) and has a higher viscos¬ 
ity (a = 0.05), we would anticipate that a strong external torque 
would be required to tear the disc, and we observe their disc does 
tear at a smaller radius than any of ours. 

4.2.3 Location of tearing radius 

The disc is expected to tear when the Lense-Thirring torque is 
larger than the internal communication in the disc. If the internal 
communication in the disc is governed by viscosity, the torques 
given in Section [2~2l can be used to estimate the upper breaking 
radius given in Equation [9] However in our simulations the disc 
internal dynamics are dominated by pressure rather than viscosity, 
hence Eauationll llmav be more appropriate. 

As the Fense-Thirring torque has a radial dependence, it is 
largest in the inner most parts of the disc and it is reasonable that 
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Figure 10. Disc tearing for different spin and viscosity combinations, with the same initial tilt of 30° and 10 5 particles shown in green (wide dashes), 10® 
particles shown in red (dashes) and 10 7 particles shown in black (solid). At the lowest resolution we observe the Bardeen-Petterson effect complete with a 
smooth transition for most discs. Increasing the resolution results in disc tearing (Bardeen-Petterson alignment) independent of our choice of viscosity or spin. 
As we do not include the affect of Einstein precession, we do not observe radial tilt oscillations in these discs. 


these discs will break at a radius smaller than predicted by Equa¬ 
tion^ Figure[l3lshows a comparison of the estimated break radius 
for the simulations inclined at 30° compared to the prediction from 
Equation [9] (upper line; assuming that a = 0.02, the average for 
our simulations) and from Equation QT| (lower line). We find that 
the disc does break at radii lower than our prediction from the vis¬ 
cous torques alone, and that the breaking radius is intermediate be¬ 
tween the predictions from Equations [9] and QT| indicating that the 
torques in our discs lie between these two extremes. The increasing 
uncertainties at low spin correspond to the decreasing convergence 
of our simulations due to mass accretion at the inner edge, seen in 
Figure [TO] 


4.2.4 Width of the rings 

The rings that are torn off during the simulations appear to be 
muc h wider than those found in the diffusive regime jNixon et alj 
l2012bl) . some up to A R/H(R) ~ 25 (where AR represents the 
ring width). It is possible for rings to form when the disc is able 
to break and differential precession is present, such as when the 
disc is subjected to Lense-Thirring precession. We therefore expect 
the width of the ring to be determined by a relative comparison be¬ 
tween the sound crossing and precessional timescales in the disc. 
We can approximate this by letting A R be the distance that a wave 
can travel in a precession time such that 


The discrepancy between the predicted and the observed 
breaking radius appears to occur at all inclinations. Using Equa- 
tion[9] the breaking radius for the 60° disc is found to be Rbreak ~ 
41f?i n which is greater than -Rout - However this disc is observed to 
break (at R < 18Ri n ), in line with the results of Figure [l3l If we 
now consider the 15° disc, it is predicted to break at Rbreak ~ 
18Rin but from the simulation we do not observe tearing. This 
could occur if the actual tearing radius is less than R; n , consistent 
with the previous results. 


/ 


— dR oc t p , 

Cs 


(32) 


across the ring. If we assume that the inner edge of the ring is at 
Rin = R — AR/2, the outer edge at R ou t = R + AR/2 and use 
the expression for the sound speed, we get 
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Figure 11. 3D renderings of discs that were initially misaligned with the black hole spin at various angles, with each simulation using 10 7 particles and shown 
after ~ 1500 orbits. The inability of the discs inclined by more than 6 > 45° to communicate the Lense-Thirring precession causes the formation of discrete 
rings which ‘tear’ and precess effectively independently before undergoing direct cancellation of angular momentum and rapid accretion. The black hole spin 
in each of these images is vertical with respect to the page (i.e. along the 2 axis). The same density scale is used as in Figure[8] 
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Figure 12. Instantaneous mass accretion of the 0°, 15°, 30°, 60°, 90°, 
120° and 150° discs run with 10 6 particles and time measured in orbits. 
The bin width is 10 times the orbit timescale and a logarithmic scale is used 
for convenience. In line with previous results, inclining the disc to the black 
hole results in mass accreting faster by almost an order of magnitude. 



Spin (a) 


Figure 13. Comparison of the breaking radius measured from the discs in¬ 
clined at 30° with our prediction of -Rbreak ( u PP er ) found by considering 
the torques in the disc (Equation [9j and inbreak t (lower) by comparing the 
sound crossing and precession timescales (Equation I Ilk 

where R is the radius that a ring of thickness A R occurs at. Fig¬ 
ure M compares the width of the rings measured from the simu¬ 
lations to this prediction. Although there are large uncertainties in 
the measurements the general trend of increasing ring width with 
R is reproduced. 

4.3 Can the disc accrete misaligned? 

Previous simulations of tilted accretion discs in the wavelike 
regime have not identified disc tearing when the disc is subjected 
to Lense-Thirring precession. The results of these thicker discs 
have found that the disc warps in the inner region, with a non- 
zero tilt at the inner edge, and then prece sses as a solid body 
dFragile & Anninosll2005luFragile et all2007h . To examine this be¬ 
haviour, we conduct a single simulation of a thick disc. We use the 
same parameters as PS9, but with an aspect ratio four times the 



Figure 14. The solid line shows the expected width of each ring of gas tom 
off in our simulations, calculated by comparing the precession timescale to 
the distance that the wave can travel (Hquation 1 33 1 . The circles show the 
measured ring widths from the simulations, where black circles indicate 
short lived rings and red circles rings that are stable for more than ~ 20 
orbits. 


initial value, such that H/R. — 0.2 a t the inner edge. This disc is 
similar to the simulation of lFragile et al.1 i2007t ). except that it has 
twice the initial tilt (and does not include magnetic fields). 

The timescales for this disc are shown in the right of Fig¬ 
ure ED Comparison with the timescales from PS9 (left) shows 
that although the precession timescale has not changed appreciably, 
the viscous and sound crossing timescales have decreased substan¬ 
tially. In the outer half of the disc we note that the viscous timescale 
is the shortest, allowing the material located there to be accreted to 
the inner regions faster than it can align. This leads to material be¬ 
ing accreted before it can align with the spin of the black hole, 
causing a non-zero tilt at the inner edge of the disc. A comparison 
between the mass accretion of this disc and our thinner PS9 simu¬ 
lation shows that there is more mass accreted by the thicker disc. 

Sim ulati ng to approximately the same time as quoted by 
iFragile et alj l2007f ). we observe the thick disc to warp in the in¬ 
ner regions but not to tear. At this time in our thin disc simulations 
we also do not observe tearing, so we continue the simulation un¬ 
til approx imately 200 orbits according to the time units specified 
bv lFragile et alj j2007h (10 times longer than their lower resolution 
simulation). The results at this time are shown in Figure[T6landl 171 
We do not observe the large i ncrease in the disc tilt at the inner 
edge that was found bv lFragile et alj j2007h (see their Figure 12), 
however in their paper this is attributed to plunging streams which 
we also do not observe. Presumably this is due to our use of the 
post-Newtonian approximation in Equation[25]and consistent with 
the discussion in Section [3.1.3l 

As the disc is four times thicker than our simulation PS9, v 
increases by a factor of 16 (even though a does not change). This 
increases the internal torque in the disc by the same factor (see 
Equation |7J, but the external torque applied is the same as for our 
disc. This should make it much harder to tear the disc, and when we 
calculate the breaking radius using Equation[9]we find that it would 
be R/Ri n ~ 3, inside the region where misaligned accretion is 
occurring. Indeed, from our results in Section l4.2.4l we would not 
expect this disc to tear at all. 
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Figure 15. The timescales in our PS9 simulation (left) and a disc that is four times thicker (right). The precession timescale does not change much between 
the thin and thick discs, however the sound crossing and viscous timescales decrease as the disc becomes thicker. The decrease in viscous time means that the 
thicker disc is able to accrete misaligned material, preventing the development of a steady state. 


5 DISCUSSION 


Despite using up to 10 7 particles, the simulations that have been 
presented are not yet c onverged. As shown in Figur e [Hand compar¬ 
ison of our results with lNelson & Papaloizoul d2000l ). increasing the 
resolution strongly affects the behaviour the disc displays. How¬ 
ever, features like disc tearing and radial oscillations are present 
in both the medium and high resolution simulations and so we 
can draw conclusions about the qualitative behaviour. Additionally, 
whilst increasing the resolution decreases the breaking radius, it 
does so by a smaller amount each time, so we are confident that the 
measured tearing radii for our non-linear simulations is an upper 
limit and that our results are close to being converged. 


The main poi nt of discussion is why o ur results differ 
to those found by iNelson & Papaloizoul d2000h . The two main 
factors are the numerical resolution and the viscosity parame- 
ter 0av- Simulations we p erformed at comparable resolution to 
[Nelson & Papaloizoul d200u) showed similar behaviour — namely 
a smooth transition between the aligned and misaligned re¬ 
gions. However, when we increased the resolution we found 
that the behaviour changes and these discs tear into two discon¬ 
nected sections (the main criterion being to adequately resolve 
the disc scale height) . This implies that low resolution prevented 
INelson & Papaloizoul J2000l) from observing disc tearing. How¬ 
ever, we also showed that the inclusion of a p viscosity, even at low 
resolution, recovers steady-state oscillations in the tilt of the disc 
midplane with respect to the blac k hole spin axis simi l ar to those 
predicte d by the linear theory of llvanov & Illario nov d 1 9 971) and 
[LOP02[ This is in contrast to the findings in Nelson & Papaloizot] 
d2000l) . where it was suggested that the tilt oscillations were short 
wavelength features which could be damped out by non-linear ef- 
fects. A s shown by our 15° simulation and in agreement with 
ILOP02L we found that the wavelength of the radial oscillation is 
of the order of the radius and is not damped out by such effects. 


The tilt oscillations that were found at linear inclinations do 
not match to the description of the ID code bv lLOPOl. and increas¬ 
ing the resolution does not reduce the discrepancy. The difference 
is likely due to the ID code assuming that the viscous timescale is 
negligibly large. In Section l4~3l it is found that the mass accretion is 
not necessarily negligible, as for discs with a larger aspect ratio we 
found it is possible for the material to accrete to the inner regions 


of the disc faster than it is able to align. This causes the disc to ac¬ 
crete misaligned material, which prevents a steady state from being 
formed and confirms that it is not possible to produce a tilt profile 
such as that described by the Bardeen-Petterson effect if the viscous 
time is too short (as predicted bv lLodato & Pringlel20061 ). Recently 
the thinnest discs in relativisti c simulations have been completed by 
iMorales Teixeira et aljd2014l ). with H/R = 0.08. Their retrograde 
simulation showed partial alignment at the inner edge, but their pro¬ 
grade simulations displayed an inner edge tilt that was greater than 
the initial condition. It is also noted that the strength of the tilt oscil¬ 
lations depends on the disc thickness, and so thick discs (and tori) 
would display weak oscillations. 

Perhaps the main caveat of our simulations is that we use an a 
viscosity to model the discs. Whilst a comparison between a purely 
hydrodynamical disc (with no explicit viscosity) and one where the 
viscosity is controlled by the MRI has shown that the behaviour 
of the disc is largel y controlled by the hydrodynamic evolution 
( ISorathia et alj|2013f) . for a complete picture of the disc evolution 
we should include magnetic fields to self-consistently generate a 
turbulent viscosity through the MRI. However it is not yet clear 
how the MRI will respond in the presence of a warp, especially at 
large angles. Additionally, we assume that our discs are vertically 
isothermal. Heat ing o f the disc due to warping may further compli¬ 
cate this picture dOgilviell2003h . 

For tilt oscillati ons and efficient wave t r ansport to occur we 
requi re H/R > a dPapaloizou & Linll 19951 : llvanov & Illarionovl 
1 19971) . Black hole accretion discs are often expected to be ge¬ 
ometrically thin and have a ~ 0.1 dKing et all [200 7). How¬ 
ever, for discs which are accreting either at very sub-Eddington 
(< 0.1Z/Edd) or near-Eddi ngton (> I/Edd) rates, the disc may be¬ 
come geometrically thick dNaravan & Quataertll2005t) . So the sim¬ 
ulations presented here may be most relevant to the low luminos¬ 
ity state of X-ray binaries where the disc can be thick and a may 
b e signific a ntly smaller t han its u sual outburst value dSm al l 19841 : 
iMever & Mever-HofmeisteJ 19841) . They are also relevant to AGN 
accreting at rates greater than Eddington and to the discs formed 
in tidal disruption events where the initial star orbit can be highly 
misaligned and the disrupted material infalling at super-Eddington 
rates. 
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Figure 16. Our thick disc simulation, similar to that of lFragile etal] j2007l) 
except that it is initialised at 30° and run for ten times longer. This disc 
is not observed to tear, as expected, but warping is observed in the inner 
regions and higher mass accretion than our thin disc. This figure is shown 
with the same density scale as Figureslsland ll II 



Figure 17. The final tilt profile of our thick disc simulation. Misaligned 
accretion occurs at the inner edge, causing a non-zero tilt and preventing a 
steady state from being formed. The inner edge features are not steady like 
the results of Figure|5] 


6 CONCLUSION 

In this work we have re-examined the Bardeen-Petterson effect in 
3D using hydrodynamical simulations of accretion discs subject to 
Lense-Thirring precession, in the regime where warps propagate in 
a wavelike manner (a < H/R). Our detailed conclusions are as 
follows: 

(i) The Bardeen-Petterson picture of an aligned inner disc 
smoothly connected to a misaligned outer disc occurs only at low 
inclinations and only when Einstein precession is not accounted 
for. Using high resolution calculations, we find both steady state 
oscillations in the disc tilt (when Einstein precession is included) 
and that discs break when they are relatively thin and highly mis¬ 
aligned to the black hole spin. 

(ii) We recover steady tilt oscillationsjbr the first time in a 3D 
hydrodynamics co de, as p redicted bv lLOPO^ . However, as the ID 
code developed bv lLQP02l assumes that mass accretion is negligi¬ 
ble, discrepancies remain between the predicted tilt profile and our 
3D results. 

(iii) Tilt oscillations are also present at higher inclinations (15°), 


showing that non-linear effects do not necessarily damp this be¬ 
haviour. 

(iv) Disc ‘tearing’ or ‘breaking’, rather than a smooth transition 
between spin-aligned and spin-misaligned parts of the disc, appears 
to be an inevitable outcome for accretion discs inclined to the black 
hole spin by more than a few degrees. This occurs regardless of 
whether the propagation of bending waves is governed by pressure 
forces or viscous stresses. 

(v) Tearing of the disc leads to rings that precess effectively 
independently. As in the diffusive regime, this can lead to direct 
cancellation of angular momentum and hence faster accretion. The 
main difference in the wavelike regime is that the rings are wider, 
with the width determined by the ratio of precession to sound cross¬ 
ing time rather than the disc scale-height. 

(vi) The Bardeen-Petterson effect cannot occur in discs where 
the viscous time is comparable to the alignment time. In this case 
the disc material is accreted misaligned. Hence it is possible to have 
discs that are misaligned with respect to the black hole spin even 
in the absence of tilt oscillations, but this can only occur at high M 
(i.e. for thick discs). 

(vii) Mass accretion rates can be enhanced by an order of mag¬ 
nitude or more when the disc is inclined with respect to the black 
hole spin. This occurs regardless of whether the disc is thick or thin. 
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APPENDIX A: MEASURING THE DISC PRECESSION 

Here we outline how the precession in the disc was measured from 
the simulations and show some example results. In order to analyse 
the properties of the discs from the simulations, we discretise the 
disc into a set of thin spherical annuli and average the properties 
of interest across the particles in each of these bins. Thi s proc ess is 
described in detail in Section 3.2.6 of iLodato & Price! d2010h . The 
twist, 'y(R), in our disc at a given radius is found by considering 
the unit angular momentum vectors at each radius bin in the disc. 
With l x (R), l y (R) and l z {R) being the unit vectors in the Cartesian 
coordinate system, we assign 

(g|), (AD 

for each radial annulus. This is repeated at every time step, so that 
we have a description of the twist as a function of time at each 
radial bin in our simulations. An example of the twist in this format 
is in Figure IaTI Here the twist is increasing in the disc when the 
gradient is positive. As the disc twists through a full 27r radians, 
the twist then jumps back to zero because Equation I A II does not 
take into account the cumulative twist angle. 
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Figure Al. An example of the description of twist in the disc as a function 
of time at a given radius. Here the time is plotted in orbits measured at the 
inner edge. 


The precession time can be measured from Figure I aTI directly 
by recording how long it takes for the disc to twist all the way 
around, equivalent to finding when the twist drops back to zero. 
This can be approximated by working out the gradient of the twist 
as a function of time and then using it to calculate the precession 
time in the disc. This is equivalent to calculating 



Because this calculation has been done at each radial annulus, we 
now have the precession time as a function of the radius in the disc, 
averaged over the length of the simulation. An example of this was 
shown in Figure[3] 

The above analysis does not take the inclination of the disc 
into account. This disc was inclined at 30°, but this angle did not 
come into our expression for t v or explicitlyjn our analys is fro m 
the twist. As outlined in the derivation by lLarwood et al J d 1996h . 
the Lense-Thirring precession is independent of the inclination be¬ 
tween the disc and the black hole spin. Repeating the above analysis 
with discs at other angles confirms this. 
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